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Abstract. The harmonic formulation of Einstein's field equations is considered, 
where the gauge conditions are introduced as dynamical constraints. The differ- 
ence between the fully constrained approach (used in analytical approximations) 
and the free evolution one (used in most numerical approximations) is pointed 
out. As a generalization, quasi-stationary gauge conditions are also discussed, 
including numerical experiments with the gauge-waves testbed. The complemen- 
tary 3+1 approach is also considered, where constraints are related instead with 
energy and momentum first integrals and the gauge must be provided separately. 
The relationship between the two formalisms is discussed in a more general frame- 
work (Z4 formalism). Different strategies in black hole simulations follow when 
introducing singularity avoidance as a requirement. More flexible quasi-stationary 
gauge conditions are proposed in this context, which can be seen as generalizations 
of the current 'freezing shift' prescriptions. 



1. Harmonic formulation 



Einstein's field equations 




are usually expressed as a system of partial diff'erential equations for the space-time 
metric gab- Soon after Einstein's 1915 paper, their mathematical structure was closely 
investigated, leading to a very convenient formulation (De Bonder 1921, 1927), namely 



2 g^'^d^cd 9ab+d(aHb) = T,abH^ + 2 g'^g^^ld.gac dfgM~Tace Tbdf]-^ ^ {Tab- 2 9ab) , (2) 



where indices inside round brackets are symmetrized and we have noted 



H^^^g'^r^.^i/^dbiVgg'")- (3) 



One can now take advantage of the general covariance of the theory. Let us define 
spacetime coordinates by a set of four independent harmonic functions, namely 




where the box stands for the general-covariant wave operator acting on functions. In 
this harmonic coordinate system, the field equations ([2]) get the simpler form 

T 

n ffah = • • • -16 n {Tab - gab) , (5) 
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where the dots stand for terms quadratic in the metric first derivatives. If we look 
at the principal part (the second derivatives terms), we see just a set of independent 
wave equations, one for every metric component. The coupling comes only through 
the right-hand-side terms. 

System ([Sj is very convenient in analytical approximations. Let us assume for 
instance that the metric admits a development of the form 

9ab = Vab + h[\^ + h^^^ +■■■ , (6) 

where ft,*^"^ is the nth-order perturbation. Then, one can express ([5]) in a recursive 
way: 

V^'d!, h^:+'^ = F,,{h^^\dh^^^) r,s<n, (7) 

which can by integrated just by inverting the standard (flat-space) wave operator. 

System ([5]), however, is not equivalent to the original field equations The 
coordinate conditions (|4]) can be interpreted as first-order constraints to be imposed 
on the solutions of the 'relaxed' system ([5]). The hard point in proving the well- 
posedness of the Cauchy problem for Einstein's equations was precisely to prove 
that the harmonic constraints ^ where actually first integrals of the relaxed system 
(Choquet-Bruhat 1952). In the analytical perturbation framework, this translates 
into the fact that fulfilling the harmonic constraints at the nth-level does not imply 
the same thing at the next level. Obtaining a true solution of the Einstein equations 
implies adjusting the integration constants in such a way that 

= , (8) 

and this must be done at every order in the perturbation development. 



1.1. Numerical Relativity applications 

The relaxed system ([5]) is also very useful in numerical approximations. The usual 
practice is using explicit time-discretization algorithms. This means that the metric 
coefficients are computed at a given time slice, assuming that one knows their values 
at the previous ones. But again fulfilling the harmonic constraints ([4]) is not granted. 
Moreover, in numerical approximations one has no adjustable integration constants. 
This means that numerical errors make the contracted Christoffel symbols obtained 
from the relaxed system to depart from their assumed harmonic (zero) value: 

{relaxed)j^a _^ q ^g^ 

In this 'free evolution' approach, one can use the non-zero values © to monitor 
the quality of the simulation. This can be done by introducing a 'zero' four-vector 
as the difference between the relaxed and the harmonic (zero) contracted Christoffel 
symbols, namely 

{relaxed)j^a (harmonic) pa 'ZZ^ (-10) 

so that true Einstein's solutions would correspond to = 0, which amounts to 
fulfilling the harmonic constraints. On the contrary, allowing for PH)) . solutions of the 
relaxed system would verify a generalized version of Q , in which ([3]) must be replaced 

by 

i?" = -g'"^ - ■ (11) 
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The vector Z" provides a new degree of freedom which arises quite naturally in 
numerical simulations. A fully covariant description is provided by the 'Z4 system' 
(Bona et al. 2003) 

Rab + Za,b + Zb, a - 8^ (T„fc - T/2 g^b) , (12) 

where the semicolon stands for the covariant derivative. We can easily see that the true 
Einstein's solutions are recovered when Z° = 0. Moreover, allowing for the contracted 
Bianchi identities, the stress-energy tensor conservation implies 

g'"Z% + R''bZ^ ^0 . (13) 

It is clear that the vanishing of Z", leading to true Einstein's solutions, is a first 
integral of the 'subsidiary system' (fT3|) . which follows from the field equations just 
assuming the stress-energy tensor conservation. 



1.2. Testing coordinate conditions 

We have seen how harmonic coordinates can be preferred for the sake of simplifying the 
(approximate) integration of the field equations. But one would like to use instead the 
General Relativity coordinate freedom for simplifying the physical interpretation of 
the results. On the contrary, choosing a convoluted coordinate system can complicate 
even the simplest physical situations. 

A good example of these unphysical gauge complications is the 'gauge waves' 
testbed (Alcubierre et al 2004). The Minkowsy (flat) metric can be written in some 
non-trivial harmonic coordinate system as: 

= F{x - t){-dt^ + dx^) + dy^ + , (14) 

where F is an arbitrary function of its argument. One could naively interpret this as 
the propagation of an arbitrary wave profile with unit speed. But it is a pure gauge 
effect, because p4|) is nothing but the Minkowsky metric in disguise. A more natural 
coordinate system should be adapted to the fact that fiat spacetime is stationary, and 
this is not granted by the harmonic condition, as (jl4p shows dramatically. 

The problem of finding 'quasi-stationary coordinates' (as stationary as possible) 
in a generic spacetime has been addressed recently (Bona et al 2005a). The idea is to 
find 'almost-Killing' vector fields by means of a standard variational principle 

(JS* = , S = J L ^ d^x , (15) 

where the Lagrangian density L is given by 

L = ^ia;b)&-'''^-^(.e.,J' (16) 

(k being an arbitrary parameter), and the variations of the vector field ^ are considered 
in a fixed spacetime. The resulting Euler-Lagrange equations get the form 

[r^'' + e'^"-fcr;c5"''];fc = o. (17) 

('almost-Killing' equation), or the equivalent one 

g'^L-M + i?ah + (1 - fc) da J - . (18) 

We will consider here the particular 'harmonic' choice k = 1. Note that, in this 
case, the subsidiary system (IT51) is nothing but condition (fTSl) for the Z° vector. We 
can then interpret that the combination . b) in the Z4 system (fT2|) gets minimized, so 



Gauge and constraint degrees of freedom 



4 




Figure 1. Gauge waves simulation with periodic boundary conditions and 
sinusoidal initial data for the metric component Qxx- The left panel corresponds 
to the harmonic coordinates case. After 100 round trips, the evolved profile (cross 
marks) nearly overlaps the initial one (continuous line), which corresponds also 
with the exact solution in the harmonic case. The right panel corresponds to the 
same simulation for the quasi-stationary coordinates defined by l|19| l . We show the 
initial data and the evolved profiles after 1, 2, 5 and 100 round trips (continuous 
line). The solution clearly approaches the Minkowsky value {gxx = 1), although 
a small residual profile remains. 



that one gets as close as possible to the original Einstein system. The name 'harmonic' 
for the k = 1 choice can be justified if we take the integral curves of ^ to be the time 
lines of our coordinate system. In this 'adapted coordinate system', condition (|17p 
reads simply 

g''dtT%, = , (19) 

which is a close generalization of the harmonic coordinates condition (see Bona et al 
2005a for more details). 

In order to test the behavior of these quasi-stationary conditions, we will consider 
the 'gauge waves' form (fT4|) of the flat metric, with the following profile: 

F =1- A Sin{ 2t:{x ~ t) ) A = 0.1, (20) 

so that the resulting metric is periodic and we can identify for instance the points 
—0.5 and 0.5 on the x axis. This allows to set up periodic boundary conditions in 
numerical simulations, so that the initial profile keeps turning around along the x 
direction. Note that the coordinate conditions (fT7|) require the use of some damping 
term in order to ensure the stability of the solutions (see Bona et al 2005b for details). 

The results of the numerical simulations are displayed in Fig. [TJ The left panel 
shows the harmonic case, where the numerical results follow very closely the exact 
solution. Only a small amount of numerical dissipation is visible after 100 round 
trips (we are using here a third-order-accurate finite volume method in order to get 
rid of the dominant dispersion error). The right panel shows the behavior of the 
quasi-stationary condition p9|) for the same initial data. We see that the amplitude 
is quickly decreasing, so that we get very close to the stationary (Minkowsky) value 
Qxx = 1 after only 5 round trips. Although a small residual profile remains, even after 
100 round trips (continuous line), the initial amplitude has been reduced by an order 
of magnitude. 
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2. 3+1 formulation 

About thirty years after the advent of the harmonic formalism, the 'evolution 
formalism' (Lichnerowicz 1944, Choquet 1956 provided a breakthrough in the 
understanding of Einstein's equations structure. Four-dimensional spacetime was 
described as sliced by constant-time hypersurfaces. Space-time dynamics was then 
described as the time evolution of the three-dimensional geometry of these surfaces. 
This '3-1-1 formalism' was very convenient at that time for studying the initial value 
problem and it is widely used today in numerical applications. 

The four-dimensional metric can be adapted to the 3-t-l geometry in the following 

way: 

ds^ = gab dx^dx^ = -o? dt^ + 7y {dx' + 13' dt) {dx^ + dt), (21) 

where 7^ is the three-dimensional metric on every slice. The 'lapse function' a 
measures the proper-time versus coordinate-time ratio when moving along the lines 
normal to the slices (/3* — 0). We will see how to take advantage of this degree of 
freedom in what follows. The 'shift' /3* measures in turn the deviation between these 
normal lines and the time lines. It does not affect the geometry of the slicing in any 
way. The extrinsic curvature of the slices is the proper time derivative of the space 
metric along the normal lines, namely: 

K,,=-\/2a{dt-Cp)-f,,, (22) 

where C stands for the Lie derivative. 

In this framework, the set of ten Einstein's field equations ([T|) can be decomposed 
into two different subsets (see Choquet 1967 for details). The space components 
provide six first-order evolution equations for Kij (amounting to second-order 
evolution equations for 7y ). The components 

rifc (G"^'' = 8 TT T"^) (23) 

(where Gab = Rab ~ 5 ^ 9ab is the Einstein tensor and Ua is the field of unit normals to 
the time slicing) provide instead four constraint equations for the pair 7^ , Kij . These 
'energy-momentum constraints' arise just from the field equations, independently of 
the coordinate gauge. Moreover, one does not get here any evolution equation for 
the lapse or the shift, which are just kinematical quantities. In order to determine 
them, one must specify four coordinate conditions, which must be added to the field 
equations in order to complete the evolution system. One can provide for instance 
four extra evolution equations 

{dt - I3^du)a = -a'Q , {dt - /3''dt)P' = -aQ' , (24) 

where Q and can be freely specified. 

In order to see the relationship between the harmonic and the 3-1-1 formalisms, 
let us decompose the contracted Christoffel symbols = g^'^T'^^, namely: 

HaT'' = aT° = Q -trK , a T, = Q, - d^a + a '■^^T, . (25) 

It follows that fixing the value of F" amounts to provide the evolution equation for the 
lapse (which determines the time slicing) , whereas fixing the value of F^ amounts to 
provide the evolution equation for the shift (which determines the time lines for a given 
slicing). The main difference is that the coordinate conditions ^ where introduced 
as constraints in the harmonic formalism, whereas the corresponding 3-1-1 conditions 
(|24p are part of the evolution system. This will have important consequences in what 
follows. 
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2.1. Gravitational collapse scenarios 

Gravitational collapse in General Relativity can lead to the arising of spacetime 
singularities, even from regular initial data (a massive star for instance). This is 
a serious issue in Numerical Relativity, as the computation can not proceed beyond 
singularity formation, even if the outside spacetime regions are regular. The lapse 
degree of freedom can be used in these cases to allow the computation to continue 
(in coordinate time) by locally diminishing the lapse (and then the proper time flow) 
in the collapsing regions, where the space volume element ^/j is diminishing. When 
done properly, the time slices do not reach the collapse singularity in a finite amount 
of coordinate time (singularity avoidance) . 

Let us illustrate this behavior by considering the harmonic slicing condition. 
Allowing for (P5|) . this means to take Q — trK. As far as singularity avoidance is 
a geometrical property of the slicing, the value of the shift is irrelevant, so we will use 
normal coordinates (zero shift) for simplicity. This condition can be easily integrated: 
allowing for the definitions ((22l [24)1 we get 

dt iVl/a) = . (26) 

It follows that the lapse tends to vanish locally (collapse of the lapse) where the space 
volume element goes to zero. Although one gets in this way arbitrarily close to the 
singularity, it can be shown that one does not actually reach it in a finite amount of 
coordinate time (Bona and Masso 1988). 

The singularity avoidance of the harmonic slicing is really a limit case, vulnerable 
to numerical errors. This is a problem in gravitational collapse scenarios, specially 
in the harmonic formalism, where we must remember that the harmonic conditions 
where not part of the relaxed system: they were just dynamical constraints. This is 
why numerical codes based in the harmonic formalism currently excise the singularity- 
forming regions out of the computational domain. This amounts to set up an artificial 
internal boundary in the strong field region. Moreover, this must be a dynamical 
boundary, capable of changing shape or moving across the numerical grid. This is a 
difficult issue, but still feasible. This has been achieved for instance in the most recent 
Numerical Relativity breakthrough, where a binary black hole simulation lasted long 
enough for extracting a consistent gravitational wave signal for the first time (Pretorius 
2005). 

An obvious alternative to dynamical excision is to consider time slicing conditions 
with stronger singularity-avoidance properties. This can be done by generalizing the 
harmonic slicing condition in the following way (Bona et al 1995): 

Q = ftrK, (27) 

where / is an arbitrary function of the lapse, so that the original harmonic slicing is 
recovered for / = 1. Most of the current binary- black-hole simulations in the BSSN 
formalism actually use the slicing condition (|27|) with the choice / — 2 /a. This 
corresponds to the 'l-|-log' slicing condition (Bernstein 1993), the name coming from 
the resulting form of the lapse in normal coordinates (zero shift): 

a = ao + ln{j/jo) ■ (28) 
It follows from that the coordinate time evolution stops before even getting close 
to the collapse singularity. It is easy to see in this case that the time evolution has 
a fixed point where the lapse vanishes (a limit surface). In normal coordinates, this 
happens when 

7 = 70 exp{-ao) , (29) 
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that is well before the vanishing of the space volume element (the initial lapse value 
being usually close to one). The appearance of a limit surface provides a safety margin 
for black hole simulations, as far as numerical errors have no chance to result into 
hitting the singularity. 

Unfortunately, singularity-avoidance does not come for free. Time slices get 
distorted in the process, mainly by a length increase along the radial direction (slice 
stretching, see Reimann and Briigmann 2004), which is nevertheless compatible with 
the overall collapse. This causes a progressive loss of resolution that can produce high- 
frequency noise in numerical simulations where the grid gets too coarse to resolve 
the profiles of the dynamical fields. This is analogous to the well known Gibbs 
phenomenon, and usually leads to code crashing as far as the noise starts growing 
quickly. 

Slice stretching is inherent to singularity avoidant slicing conditions, but the 
appearance of high-frequency noise can be delayed, even avoided, in many ways. 
One can just increase the grid resolution or use instead more specific tools: artificial 
dissipation terms (Gustavson et al 1995) or advanced finite volume methods, like the 
ones currently used in Computational Fluid Dynamics (for a recent implementation, 
see Alic et al 2007). Just after Pretorius paper, many groups published improved 
binary-black-hole simulations, obtained with 3-1-1 codes using the 'l-|-log' slicing 
condition (Campanelli et al 2006, Baker et al 2006, Diener et al 2006). 



2.2. Quasi- stationary gauge conditions and singularity avoidance 

At this point, it is interesting to ask whether the quasi-stationary coordinate conditions 
derived from the variational principle (|15p can be made compatible with the singularity 
avoidance requirement. Note that, when adapting our spacetime coordinates to a 
solution <^ of the almost-Killing equation ^T7\ . we are demanding two different things: 
the time lines are chosen to be the integral curves of ^ and the time coordinate is 
chosen to be the preferred afhne parameter associated to these lines. Although the first 
requirement fits perfectly into the idea of getting a quasi-stationary gauge condition, 
the second one lacks of a clear physical motivation. It seems that singularity avoidance 
is not enforced in this way. 

A better strategy is to choose a priori the time coordinate, that is a spacetime 
slicing given by 

(/)(x°) — constant , (30) 

with a view to enforce singularity avoidance. Then, we can use this time coordinate 
as a parameter along the integral lines of the almost-Killing vector ^ , by requiring 

r - 1 . (31) 

This amounts to constraint the vector ^ to fulfill (|3ip in the minimization process. 
In other words, we introduce a Lagrange multiplier and minimize 

L' = L + X (r9„ <^ - 1) (32) 

(for any (jj given a priori) instead of the original Lagrangian (|16|) . 

The resulting Euler-Lagrange equations include now the constraint (I3ip . plus the 
system 

[r-' + i'"''-ke;c9"'];b = Xd^<P . (33) 
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Figure 2. Plots of the lapse function evolution for two splierically symmetric 
(ID) black-hole simulations. The 1+log slicing condition has been used in both 
cases, producing the lapse collapse. The space coordinates, however, are different: 
we compare normal coordinates (no shift) with quasistationary coordinates (AKE 
shift). The shift delays the collapse, as it provides an outgoing speed for the grid 
nodes. A smoothing of the slopes can also be seen as a side effect. The logarithmic 
character of the grid makes these effects less apparent for larger values of the r] 
coordinate. 

which generalizes the almost-KiUing equation (fT7|) . Using adapted coordinates, 

= i , ^ = ft , (34) 
it can be written as 

g'^dtV^^ + (1 - fc) g'^'dtT'^bc = A 5'^* . (35) 

Let us now split condition (|35p into its 3+1 components. The time component 
can be ignored, because it just provides the value of the Lagrange multiplier itself. 
Remember that the time slicing was chosen a priori, so that we do not need any further 
condition for the lapse. The (downstairs) space component, 

g.ag'''dtT%^ + il-k)dtT^,, = 0, (36) 

provides a (second order) evolution equation for the shift. In this way, we have 
completely splitted the gauge conditions: the shift equation ([36| should work with 
any time slicing, that can be freely chosen a priori. 

The quasi-stationary shift condition (j36p is completely independent of the value 
of the Lagrange multiplier. This means that we would get the same condition from 
the original unconstrained Lagrangian. We can conclude that the slicing constraint 
(|3ip does not affect the minimization process in the shift sector. 

Let us test these conditions in a spherically symmetric black-hole simulation. We 
will write the Schwarzschild line element in the 'wormhole' form: 

ds^ = -( tanhrj ^ dt^ + im^ ( coshrj/2 f ( drj^ + dfl^ ) , (37) 

which can be obtained from the isotropic form by the following coordinate 
transformation 

r — m/2 exp{ r/ ) . (38) 

We will combine the '1+log' slicing prescriptio with the quasi-stationary shift 
condition ([36]) with k — 1/2. We run the simulation up to lOOOm, as in the zero 
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shift case (normal coordinates). We see in Fig. [2] that the lapse collapses as usual, 
avoiding the singularity in both cases as expected. The effect of the shift is adding 
some outgoing speed to the grid nodes, so that the advance of the collapse front across 
the grid is delayed. We have added to the shift condition a standard damping 
term (Bona et al 2005b) in order to avoid the shift values to grow without control. 
Note that the slopes at the collapse front are smoothed out, so they can be better 
resolved. The logarithmic character of the grid makes these effects less apparent for 
larger values of the 77 coordinate. 

Our results prove that the quasi-stationary shift conditions ([55]) are actually 
compatible with standard singularity-avoidant slicing conditions. We have found, 
however, that the results can depend crucially of particular value of the parameter k. 
It is suggestive that, at least in the particular case shown here, the best results are 
obtained with the k — 1/2 choice. This is a very special value, because the minimum 
principle (jlSp leads in this case to a minimisation of the conformal-Killing equation: a 
quasi-conformal shift condition. This opens an interesting perspective for future work. 
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